
## load data
surv <- read.csv("data/yougov experimental data.csv")

## rename questions
names(surv)[c(
  2,
  3:8,
  9:12,
  13:15,
  16:18
  )] <- 
  c("treatment_arm",
    "thrm_cong","thrm_corp","thrm_ent","thrm_cath","thrm_who","thrm_nfl",
    "TradeOpenFavorTA1","TradeDealFavorTA1","TradeOpenFavorTA2","TradeDealFavorTA2",
    "med_job","med_inq","med_pol",
    "employment","emp_size","emp_size_unemp"
   )
head(surv)

# treatment indicator variable: 
# Treatment Arm 1 is the "big firms" treatment; Treatment Arm 2 is the "competitive industries" treatment
surv$treatment <- NA
surv$treatment[surv$treatment_arm == "Arm 1"] <- 1
surv$treatment[surv$treatment_arm == "Arm 2"] <- 0

# support for trade and trade agreements
surv$TradeOpenFavor <- NA
surv$TradeOpenFavor[surv$treatment == 1] <- surv$TradeOpenFavorTA1[surv$treatment == 1]
surv$TradeOpenFavor[surv$treatment == 0] <- surv$TradeOpenFavorTA2[surv$treatment == 0]
surv$TradeOpenFavor <- factor(surv$TradeOpenFavor, levels = c("Oppose a great deal","Oppose somewhat",
  "Neither favor nor oppose","Favor somewhat","Favor a great deal"))
surv$TradeOpenFavor <- as.numeric(surv$TradeOpenFavor) 

# support for trade and trade agreements
surv$TradeDealFavor <- NA
surv$TradeDealFavor[surv$treatment == 1] <- surv$TradeDealFavorTA1[surv$treatment == 1]
surv$TradeDealFavor[surv$treatment == 0] <- surv$TradeDealFavorTA2[surv$treatment == 0]
surv$TradeDealFavor <- factor(surv$TradeDealFavor, levels = c("Oppose a great deal","Oppose somewhat",
  "Neither favor nor oppose","Favor somewhat","Favor a great deal"))
surv$TradeDealFavor <- as.numeric(surv$TradeDealFavor) 

# employer size
surv$emp_size_comb <- surv$emp_size; 
surv$emp_size_comb <- factor(surv$emp_size_comb, levels = c("1-5","6-19","20-49","50-199",
  "200-999","1000-9999","More than 10,000","I have never been employed"))
surv$emp_size_comb[is.na(surv$emp_size_comb)] <- surv$emp_size_unemp[is.na(surv$emp_size_comb)]

# race
race2 <- rep("Other non-white", nrow(surv))
race2[surv$race == "White"] <- "White"
race2[surv$race == "Black"] <- "Black"
race2[surv$race == "Hispanic"] <- "Latino"
race2[surv$race == "Asian"] <- "AAPI"
surv$race <- factor(race2, levels = c("White","Black","Latino","AAPI","Other non-white"))

# gender
surv$gender <- factor(surv$gender, levels = c("Female","Male"))

# mediators convert to numeric.
surv$med_job <- factor(surv$med_job, levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree"))
surv$med_job <- as.numeric(surv$med_job)
surv$med_inq <- factor(surv$med_inq, levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree"))
surv$med_inq <- as.numeric(surv$med_inq)
surv$med_pol <- factor(surv$med_pol, levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree","Somewhat agree","Strongly agree"))
surv$med_pol <- as.numeric(surv$med_pol)

# party and ideology, turn to numeric scales
surv$party <- factor(surv$pid7, levels = c("Strong Democrat","Not very strong Democrat","Lean Democrat","Independent",
  "Lean Republican","Not very strong Republican","Strong Republican"))
surv$party <- as.numeric(surv$party)
surv$ideology <- factor(surv$ideo5, levels = c("Very liberal","Liberal","Moderate","Conservative","Very conservative"))
surv$ideology <- as.numeric(surv$ideology)
surv$income <- factor(surv$faminc_new, levels = c("Less than $10,000","$10,000 - $19,999","$20,000 - $29,999",
           "$30,000 - $39,999","$40,000 - $49,999","$50,000 - $59,999","$60,000 - $69,999","$70,000 - $79,999",
           "$80,000 - $99,999","$100,000 - $119,999","$120,000 - $149,999","$150,000 - $199,999","$200,000 - $249,999",
           "$250,000 - $349,999","$350,000 - $499,999","$500,000 or more","Prefer not to say"))
surv$income <- as.numeric(surv$income)
surv$income[surv$income==17] <- NA

# age
surv$age <- as.numeric(surv$birthyr)

# collapse education to numeric binary
surv$college <- as.numeric(surv$educ %in% c("Some college","2-year","4-year","Post-grad"))

# collapse employed into smaller number of categories and make numeric
surv$employed <- as.numeric(surv$employment %in% c("Full-time","Part-time"))
surv$unemployed <- as.numeric(surv$employment %in% c("Temporarily laid off","Unemployed"))
surv$retstudis <- as.numeric(surv$employment %in% c("Retired","Student","Permanently disabled","Homemaker","Other"))

surv$emp_size_comb <- factor(surv$emp_size_comb, c("1-5","6-19","20-49","50-199","200-999","1000-9999","More than 10,000"))
surv$emp_size_comb <- as.numeric(surv$emp_size_comb)

# convert empsize into a above/below median dummy
surv$emp_size_big <- as.numeric(as.numeric(surv$emp_size_comb) > median(as.numeric(surv$emp_size_comb), na.rm = TRUE))

# convert used thermometer into above/below median dummies.
surv$thrm_corp_pos <- as.numeric(as.numeric(surv$thrm_corp) > median(as.numeric(surv$thrm_corp), na.rm = TRUE))

##############################
## Average treatment affect ##
##############################

mod1 <- lm(TradeOpenFavor ~ treatment, data = surv, weights = weight);
mod2 <- lm(TradeOpenFavor ~ treatment + age + gender + race, data = surv, weights = weight); 
mod3 <- lm(TradeOpenFavor ~ treatment + age + gender + race + college + income + employed + unemployed, data = surv, weights = weight); 
mod4 <- lm(TradeOpenFavor ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight); 

# coef(mod1)[2]/sqrt(wtd.var(surv$TradeOpenFavor, weights = surv$weights))

mods <- list(summary(mod1)$coef, summary(mod2)$coef, summary(mod3)$coef, summary(mod4)$coef)
vars <- c(rownames(mods[[4]])[2:14], "(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)","Intercept") 
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(28,29); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB4.tex")

# CI's
LBs <- c(coef(mod1)[2] - qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] - qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] - qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] - qnorm(.975)*summary(mod4)$coef[2,2])
UBs <- c(coef(mod1)[2] + qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] + qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] + qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] + qnorm(.975)*summary(mod4)$coef[2,2])
CIs <- apply(cbind(LBs,UBs), 2, myround, 2)

est <- tab[1,]
ci <- paste("\\multicolumn{1}{c}{", paste("$[", CIs[,1], ",", CIs[,2], "]$\\phantom{'''''}", sep = ""), "}", sep = "")
ci <- gsub("-","\\\\sminus", ci); ci <- gsub("0\\.","\\.", ci);
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "\\phantom{''''}}", sep = ""); tab <- rbind(tab, N)
tab <- rbind(est, ci, N)
rownames(tab) <- c("Average treatment effect", "ATE 95\\% CI", "N")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(cbind(rownames(tab), tab), "paper/tables/Table1_top.tex")

# 
mod1 <- lm(TradeDealFavor ~ treatment, data = surv, weights = weight); 
mod2 <- lm(TradeDealFavor ~ treatment + age + gender + race, data = surv, weights = weight); 
mod3 <- lm(TradeDealFavor ~ treatment + age + gender + race + college + income + employed + unemployed, data = surv, weights = weight); 
mod4 <- lm(TradeDealFavor ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight); 

mods <- list(summary(mod1)$coef, summary(mod2)$coef, summary(mod3)$coef, summary(mod4)$coef)
vars <- c(rownames(mods[[4]])[2:14], "(Intercept)")
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)","Intercept") 
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(28,29); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB5.tex")

# CI's
LBs <- c(coef(mod1)[2] - qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] - qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] - qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] - qnorm(.975)*summary(mod4)$coef[2,2])
UBs <- c(coef(mod1)[2] + qnorm(.975)*summary(mod1)$coef[2,2],
  coef(mod2)[2] + qnorm(.975)*summary(mod2)$coef[2,2],
  coef(mod3)[2] + qnorm(.975)*summary(mod3)$coef[2,2],
  coef(mod4)[2] + qnorm(.975)*summary(mod4)$coef[2,2])
CIs <- apply(cbind(LBs,UBs), 2, myround, 2)

est <- tab[1,]
ci <- paste("\\multicolumn{1}{c}{", paste("$[", CIs[,1], ",", CIs[,2], "]$\\phantom{'''''}", sep = ""), "}", sep = "")
ci <- gsub("-","\\\\sminus", ci); ci <- gsub("0\\.","\\.", ci);
N <- c(sum(summary(mod1)$df[1:2]),sum(summary(mod2)$df[1:2]),sum(summary(mod3)$df[1:2]),sum(summary(mod4)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "\\phantom{''''}}", sep = ""); tab <- rbind(tab, N)
tab <- rbind(est, ci, N)
rownames(tab) <- c("Average treatment effect", "ATE 95\\% CI", "N")
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(cbind(rownames(tab), tab), "paper/tables/Table1_bottom.tex")

#####################################
## Heterogeneous treatment effects ##
#####################################

## Employmer size
mod1a <- lm(TradeOpenFavor ~ treatment*emp_size_big, data = surv, weights = weight); summary(mod1a)

## Employmer size
mod1b <- lm(TradeOpenFavor ~ treatment*I(emp_size_comb-1), data = surv, weights = weight); summary(mod1b)
betas <- mvrnorm(10000, coef(mod1b), vcov(mod1b))
emp_size_seq <- seq(0, 6, by = 1/3)
control_mat <- data.frame(int = 1, treatment = 0, emp_size_comb = emp_size_seq, interaction = 0)
control_sims <- betas %*% t(control_mat)
treated_mat <- data.frame(int = 1, treatment = 1, emp_size_comb = emp_size_seq, interaction = emp_size_seq)
treated_sims <- betas %*% t(treated_mat)

ate.means <- apply(treated_sims - control_sims, 2, mean)
ate.ub <- apply(treated_sims - control_sims, 2, quantile, .975)
ate.lb <- apply(treated_sims - control_sims, 2, quantile, .025)

# plot(ate.means ~ emp_size_seq, type = "l", ylim = c(-.7,.2))
# lines(ate.ub ~ emp_size_seq, type = "l", ylim = c(-.7,.2))
# lines(ate.lb ~ emp_size_seq, type = "l", ylim = c(-.7,.2))

## positive attitude toward corporations
mod2a <- lm(TradeOpenFavor ~ treatment*thrm_corp_pos, data = surv, weights = weight); summary(mod2a)
quantile(surv$thrm_corp, seq(0,1,.1))

## Corporate thermometer
mod2b <- lm(TradeOpenFavor ~ treatment*thrm_corp, data = surv, weights = weight);
betas <- mvrnorm(10000, coef(mod2b), vcov(mod2b))
thrm_corp_seq <- seq(0, 100, by = 5)
control_mat <- data.frame(int = 1, treatment = 0, thrm_corp = thrm_corp_seq, interaction = 0)
control_sims <- betas %*% t(control_mat)
treated_mat <- data.frame(int = 1, treatment = 1, thrm_corp = thrm_corp_seq, interaction = thrm_corp_seq)
treated_sims <- betas %*% t(treated_mat)

ate.means <- apply(treated_sims - control_sims, 2, mean)
ate.ub <- apply(treated_sims - control_sims, 2, quantile, .975)
ate.lb <- apply(treated_sims - control_sims, 2, quantile, .025)

pdf("paper/figures/FigureB1_bottom.pdf", height = 4, width = 5)

par(mar = c(2.5, 2, 2, 1)) # Set the margin on all sides to 2
plot(ate.means ~ thrm_corp_seq, type = "l", xlim = c(-1,102), ylim = c(-.78,.6), lwd = 2, 
  axes = FALSE, xlab = "", ylab = "")
lines(ate.ub ~ thrm_corp_seq, type = "l", ylim = c(-.72,.4))
lines(ate.lb ~ thrm_corp_seq, type = "l", ylim = c(-.72,.4))
library(scales)

smoother <- density(as.numeric(surv$thrm_corp), kernel = c("gaussian"), width = 2)
lines(smoother$y + (-.76) ~ smoother$x)
jits <- jitter(surv$thrm_corp, factor = 1.5)
segments(x0 = jits, x1 = jits, y0 = -.8, y1 = -.78, col = alpha("black",.2))

mtext(side=2, "Conditional ATE", cex = .6, line = 1, adj = .6)
axis(2, at = c(-.6,-.3,0,.3,.6), labels = c("-.6  ","-.3  "," 0 ",".3 ",".6 "), 
  cex.axis = .5, tck = -.01, padj = 2.75, lwd = .5) 

mtext(side=1, "Views of Corporate America", cex = .6, line = .6, adj = .5)
axis(1, at = c(0,25,50,75,100), labels = c("Negative","","Neutral","","Positive"), 
  cex.axis = .5, tck = -.01, padj = -3.5, lwd = .5) 

mtext(side=3, "YouGov Experiment: Treatment effect conditional on attitude toward corporations", cex = .7, line = .3, adj = .5)

dev.off()

pdf("paper/figures/Figure1.pdf", height = 4, width = 5)

par(mar = c(2.5, 2, 2, 1)) # Set the margin on all sides to 2
plot(ate.means ~ thrm_corp_seq, type = "l", xlim = c(-1,102), ylim = c(-.78,.6), lwd = 2, 
  axes = FALSE, xlab = "", ylab = "")
lines(ate.ub ~ thrm_corp_seq, type = "l", ylim = c(-.72,.4))
lines(ate.lb ~ thrm_corp_seq, type = "l", ylim = c(-.72,.4))
library(scales)

smoother <- density(as.numeric(surv$thrm_corp), kernel = c("gaussian"), width = 2)
lines(smoother$y + (-.76) ~ smoother$x)
jits <- jitter(surv$thrm_corp, factor = 1.5)
segments(x0 = jits, x1 = jits, y0 = -.8, y1 = -.78, col = alpha("black",.2))

mtext(side=2, "Conditional ATE", cex = .6, line = 1, adj = .6)
axis(2, at = c(-.6,-.3,0,.3,.6), labels = c("-.6  ","-.3  "," 0 ",".3 ",".6 "), 
  cex.axis = .5, tck = -.01, padj = 2.75, lwd = .5) 

mtext(side=1, "Views of Corporate America", cex = .6, line = .6, adj = .5)
axis(1, at = c(0,25,50,75,100), labels = c("Negative","","Neutral","","Positive"), 
  cex.axis = .5, tck = -.01, padj = -3.5, lwd = .5) 

mtext(side=3, "Treatment effect conditional on attitude toward corporations", cex = .7, line = .3, adj = .5)

dev.off()

# cor(surv[,c("emp_size_comb","thrm_corp")], use = "complete.obs")

mods <- list(summary(mod2a)$coef, summary(mod2b)$coef, summary(mod1a)$coef, summary(mod1b)$coef)
vars <- c(rownames(mods[[1]])[2:4], rownames(mods[[2]])[3:4], rownames(mods[[3]])[3:4], rownames(mods[[4]])[3:4])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17)] <- 
  c("Treated (Large firms ben./Small firms harmed)",
    "Positive view of corporations","Treated$\\cdot$ Pos. view of corps.","View of corporations (0-100)","Treated$\\cdot$ View of corps.",
    "Large employer","Treated$\\cdot$ Large employer","Employer size (0-6)","Treated$\\cdot$ Employer size") 
N <- c(sum(summary(mod2a)$df[1:2]),sum(summary(mod2b)$df[1:2]),sum(summary(mod1a)$df[1:2]),sum(summary(mod1b)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(18,19); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/Table2.tex")

mod1a <- lm(TradeDealFavor ~ treatment*emp_size_big, data = surv); summary(mod1a)
mod1b <- lm(TradeDealFavor ~ treatment*I(emp_size_comb-1), data = surv); summary(mod1b)
mod2a <- lm(TradeDealFavor ~ treatment*thrm_corp_pos, data = surv); summary(mod2a)
mod2b <- lm(TradeDealFavor ~ treatment*thrm_corp, data = surv);

mods <- list(summary(mod2a)$coef, summary(mod2b)$coef, summary(mod1a)$coef, summary(mod1b)$coef)
vars <- c(rownames(mods[[1]])[2:4], rownames(mods[[2]])[3:4], rownames(mods[[3]])[3:4], rownames(mods[[4]])[3:4])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17)] <- 
  c("Treated (Large firms ben./Small firms harmed)",
    "Positive view of corporations","Treated$\\cdot$ Pos. view of corps.","View of corporations (0-100)","Treated$\\cdot$ View of corps.",
    "Large employer","Treated$\\cdot$ Large employer","Employer size (0-6)","Treated$\\cdot$ Employer size") 
N <- c(sum(summary(mod2a)$df[1:2]),sum(summary(mod2b)$df[1:2]),sum(summary(mod1a)$df[1:2]),sum(summary(mod1b)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(18,19); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB7_bottom.tex")

###############
## Mediation ##
###############

## Trade openness.
set.seed(12345)
surv$treatment <- as.numeric(surv$treatment == "1")
med.fit1 <- suppressWarnings(polr(as.factor(med_job) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE, weights = weight))
out.fit1 <- lm(TradeOpenFavor ~ med_job + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight)
med.out1 <- mediate(med.fit1, out.fit1, treat = "treatment", mediator = "med_job",
  robustSE = FALSE, sims = 2000)

med.fit2 <- suppressWarnings(polr(as.factor(med_inq) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE, weights = weight))
out.fit2 <- lm(TradeOpenFavor ~ med_inq + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight)
med.out2 <- mediate(med.fit2, out.fit2, treat = "treatment", mediator = "med_inq",
  robustSE = FALSE, sims = 2000)

med.fit3 <- suppressWarnings(polr(as.factor(med_pol) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE, weights = weight))
out.fit3 <- lm(TradeOpenFavor ~ med_pol + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight)
med.out3 <- mediate(med.fit3, out.fit3, treat = "treatment", mediator = "med_pol",
  robustSE = FALSE, sims = 2000)

# 
tab3top <- rbind( 
c(summary(med.out1)$tau.coef, summary(med.out1)$tau.p, summary(med.out1)$tau.ci),

c(summary(med.fit2)$coef[1,c(1,3)], summary(med.fit2)$coef[1,1] - 1.96*summary(med.fit2)$coef[1,2], summary(med.fit2)$coef[1] + 1.96*summary(med.fit2)$coef[1,2]),
c(summary(med.out2)$d1, summary(med.out2)$d1.p, summary(med.out2)$d1.ci),
c(summary(med.out2)$z1, summary(med.out2)$z1.p, summary(med.out2)$z1.ci),

c(summary(med.fit3)$coef[1,c(1,3)], summary(med.fit3)$coef[1,1] - 1.96*summary(med.fit3)$coef[1,2], summary(med.fit3)$coef[1] + 1.96*summary(med.fit3)$coef[1,2]),
c(summary(med.out3)$d1, summary(med.out3)$d1.p, summary(med.out3)$d1.ci),
c(summary(med.out3)$z1, summary(med.out3)$z1.p, summary(med.out3)$z1.ci),

c(summary(med.fit1)$coef[1,c(1,3)], summary(med.fit1)$coef[1,1] - 1.96*summary(med.fit1)$coef[1,2], summary(med.fit1)$coef[1] + 1.96*summary(med.fit1)$coef[1,2]),
c(summary(med.out1)$d1, summary(med.out1)$d1.p, summary(med.out1)$d1.ci),
c(summary(med.out1)$z1, summary(med.out1)$z1.p, summary(med.out1)$z1.ci)
)
tab3top[,c(1,3,4)] <- apply(tab3top[,c(1,3,4)], 2, myround, 2)
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(tab3top, "paper/tables/Table3_top.tex")

sum.med.fit1 <- cbind(summary(med.fit1)$coef, 2*pnorm(-abs(summary(med.fit1)$coef[,3])))
sum.med.fit2 <- cbind(summary(med.fit2)$coef, 2*pnorm(-abs(summary(med.fit2)$coef[,3])))
sum.med.fit3 <- cbind(summary(med.fit3)$coef, 2*pnorm(-abs(summary(med.fit3)$coef[,3])))
mods <- list(sum.med.fit2, summary(out.fit2)$coef, sum.med.fit3, summary(out.fit3)$coef, sum.med.fit1, summary(out.fit1)$coef)
vars <- c(rownames(mods[[1]])[1], rownames(mods[[2]])[2], rownames(mods[[4]])[2], rownames(mods[[6]])[c(2,4:15)])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){# i <- 1
  mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))
}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Mediator: Econ. inequality","Mediator: Socio-pol. inequality","Mediator: Job insecurity",
    "Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)") 
N <- c(1511,sum(summary(out.fit2)$df[1:2]),1511,sum(summary(out.fit3)$df[1:2]),1511,sum(summary(out.fit1)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(32,33); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB11.tex")

## Trade agreements.
set.seed(12345)
surv$treatment <- as.numeric(surv$treatment == "1")
med.fit1 <- suppressWarnings(polr(as.factor(med_job) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE, weights = weight))
out.fit1 <- lm(TradeDealFavor ~ med_job + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight)
med.out1 <- mediate(med.fit1, out.fit1, treat = "treatment", mediator = "med_job",
  robustSE = FALSE, sims = 2000)

med.fit2 <- suppressWarnings(polr(as.factor(med_inq) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE, weights = weight))
out.fit2 <- lm(TradeDealFavor ~ med_inq + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight)
med.out2 <- mediate(med.fit2, out.fit2, treat = "treatment", mediator = "med_inq",
  robustSE = FALSE, sims = 2000)

med.fit3 <- suppressWarnings(polr(as.factor(med_pol) ~ treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, Hess = TRUE, weights = weight))
out.fit3 <- lm(TradeDealFavor ~ med_pol + treatment + age + gender + race + college + income + employed + unemployed + 
  party + ideology, data = surv, weights = weight)
med.out3 <- mediate(med.fit3, out.fit3, treat = "treatment", mediator = "med_pol",
  robustSE = FALSE, sims = 2000)

# 
tab3bottom <- rbind( 
c(summary(med.out1)$tau.coef, summary(med.out1)$tau.p, summary(med.out1)$tau.ci),

c(summary(med.fit2)$coef[1,c(1,3)], summary(med.fit2)$coef[1,1] - 1.96*summary(med.fit2)$coef[1,2], summary(med.fit2)$coef[1] + 1.96*summary(med.fit2)$coef[1,2]),
c(summary(med.out2)$d1, summary(med.out2)$d1.p, summary(med.out2)$d1.ci),
c(summary(med.out2)$z1, summary(med.out2)$z1.p, summary(med.out2)$z1.ci),

c(summary(med.fit3)$coef[1,c(1,3)], summary(med.fit3)$coef[1,1] - 1.96*summary(med.fit3)$coef[1,2], summary(med.fit3)$coef[1] + 1.96*summary(med.fit3)$coef[1,2]),
c(summary(med.out3)$d1, summary(med.out3)$d1.p, summary(med.out3)$d1.ci),
c(summary(med.out3)$z1, summary(med.out3)$z1.p, summary(med.out3)$z1.ci),

c(summary(med.fit1)$coef[1,c(1,3)], summary(med.fit1)$coef[1,1] - 1.96*summary(med.fit1)$coef[1,2], summary(med.fit1)$coef[1] + 1.96*summary(med.fit1)$coef[1,2]),
c(summary(med.out1)$d1, summary(med.out1)$d1.p, summary(med.out1)$d1.ci),
c(summary(med.out1)$z1, summary(med.out1)$z1.p, summary(med.out1)$z1.ci)
)
tab3bottom[,c(1,3,4)] <- apply(tab3bottom[,c(1,3,4)], 2, myround, 2)
addtorow <- list(); addtorow$pos <- list(1); addtorow$command <- c("")
ptable(tab3bottom, "paper/tables/Table3_bottom.tex")

sum.med.fit1 <- cbind(summary(med.fit1)$coef, 2*pnorm(-abs(summary(med.fit1)$coef[,3])))
sum.med.fit2 <- cbind(summary(med.fit2)$coef, 2*pnorm(-abs(summary(med.fit2)$coef[,3])))
sum.med.fit3 <- cbind(summary(med.fit3)$coef, 2*pnorm(-abs(summary(med.fit3)$coef[,3])))
mods <- list(sum.med.fit2, summary(out.fit2)$coef, sum.med.fit3, summary(out.fit3)$coef, sum.med.fit1, summary(out.fit1)$coef)
vars <- c(rownames(mods[[1]])[1], rownames(mods[[2]])[2], rownames(mods[[4]])[2], rownames(mods[[6]])[c(2,4:15)])
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){# i <- 1
  mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
  stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))
}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31)] <- 
  c("Treated (Large firms ben./Small firms harmed)","Mediator: Econ. inequality","Mediator: Socio-pol. inequality","Mediator: Job insecurity",
    "Age","Male","Black","Latino","AAPI","Other non-white","College-educated","Income","Employed","Unemployed",
    "Party (D=1,R=7)","Ideology (L=1,C=7)") 
N <- c(1511,sum(summary(out.fit2)$df[1:2]),1511,sum(summary(out.fit3)$df[1:2]),1511,sum(summary(out.fit1)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(32,33); addtorow$command <- c(paste0('\\midrule '),paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "paper/tables/TableB12.tex")



